In this script, a trimmed-down version (excluding any plots that are the same as in Script 02, such as the loadings of the cell count PCs) of Script 02 will be run for several clocks. In other words, we will test how cell counts associate with several epigenetic clocks, both DNAmAge (the clock’s age estimation) and AgeAccel (the deviation of DNAmAge from calendar age). Since we see that the PC-approach works best, we will only apply this method for the clocks.
Load sample sheet.
setwd("/exports/molepi/RSC_BIOS/Users/tjonkman/cellcounts")
.libPaths("/exports/molepi/RSC_BIOS/Users/tjonkman/Packages/4.3.1")
load("01_sample_sheet.rda")
dim(ss)
## [1] 4058 27
#Load PCA data.
load("02_pca_df.rda")
ss.pc <- cbind(ss, pca.df)
library(ggplot2)
library(reshape2)
#Print the correlations between the investigated clocks (residuals for all except dunedinPACE).
round(cor(ss[,c(5:10)]), 2)
## hannum horvath zhang phenoage grimage dunedinpace
## hannum 1.00 0.96 0.98 0.96 0.95 0.49
## horvath 0.96 1.00 0.97 0.95 0.93 0.45
## zhang 0.98 0.97 1.00 0.96 0.96 0.45
## phenoage 0.96 0.95 0.96 1.00 0.95 0.54
## grimage 0.95 0.93 0.96 0.95 1.00 0.57
## dunedinpace 0.49 0.45 0.45 0.54 0.57 1.00
round(cor(ss[,c(23:27, 10)]), 2)
## hannumRes horvathRes zhangRes phenoageRes grimageRes dunedinpace
## hannumRes 1.00 0.55 0.66 0.59 0.39 0.27
## horvathRes 0.55 1.00 0.60 0.55 0.22 0.15
## zhangRes 0.66 0.60 1.00 0.55 0.23 0.15
## phenoageRes 0.59 0.55 0.55 1.00 0.49 0.40
## grimageRes 0.39 0.22 0.23 0.49 1.00 0.60
## dunedinpace 0.27 0.15 0.15 0.40 0.60 1.00
#DNAmAge
plot.data <- ss[,c("age", "hannum", "horvath", "zhang", "phenoage", "grimage")]
colnames(plot.data) <- c("Calendar Age", "Hannum", "Horvath", "Zhang", "PhenoAge", "GrimAge")
plot.data <- melt(plot.data, id.vars = "Calendar Age", variable.name = "Clock", value.name = "DNAmAge")
library(ggpubr)
ggplot(data = plot.data, aes(x = `Calendar Age`, y = DNAmAge)) +
geom_point(shape = 1, size = 1) +
geom_abline(slope = 1, intercept = 0, linetype = "dashed", color= "red") +
geom_smooth(method = "lm") +
stat_cor(vjust = 0.5) +
facet_wrap(facets = vars(Clock), ncol = 2, scales = "free", dir = "v") +
scale_x_continuous(limits = c(18, 87), breaks = c(0, 20, 40, 60, 80)) +
scale_y_continuous(limits = c(-5, 123), breaks = c(0, 20, 40, 60, 80, 100, 120)) +
theme_bw()
## `geom_smooth()` using formula = 'y ~ x'
#AgeAccel
plot.data <- ss[,c("age", "hannumRes", "horvathRes", "zhangRes", "phenoageRes", "grimageRes")]
colnames(plot.data) <- c("Calendar Age", "Hannum", "Horvath", "Zhang", "PhenoAge", "GrimAge")
plot.data <- melt(plot.data, id.vars = "Calendar Age", variable.name = "Clock", value.name = "AgeAccel")
ggplot(data = plot.data, aes(x = `Calendar Age`, y = AgeAccel)) +
geom_point(shape = 1, size = 1) +
geom_abline(slope = 0, intercept = 0, linetype = "dashed", color= "red") +
geom_smooth(method = "lm") +
stat_cor(vjust = 0.5) +
facet_wrap(facets = vars(Clock), ncol = 2, scales = "free", dir = "v") +
scale_x_continuous(limits = c(18, 87), breaks = c(0, 20, 40, 60, 80)) +
scale_y_continuous(limits = c(-38, 68), breaks = c(-20, 0, 20, 40, 60)) +
theme_bw()
## `geom_smooth()` using formula = 'y ~ x'
#DunedinPACE
plot.data <- ss[,c("age", "dunedinpace")]
colnames(plot.data) <- c("Calendar Age", "DunedinPACE")
plot.data <- melt(plot.data, id.vars = "Calendar Age", variable.name = "Clock", value.name = "AgeAccel")
ggplot(data = plot.data, aes(x = `Calendar Age`, y = AgeAccel)) +
geom_point(shape = 1, size = 1) +
geom_abline(slope = 0, intercept = 52, linetype = "dashed", color = "red") +
geom_smooth(method = "lm", ) +
stat_cor(vjust = 0.5) +
facet_wrap(facets = vars(Clock), nrow = 1, scales = "free") +
scale_x_continuous(limits = c(18, 87), breaks = c(0, 20, 40, 60, 80)) +
scale_y_continuous(limits = c(26, 86), breaks = c(20, 40, 60, 80)) +
theme_bw()
## `geom_smooth()` using formula = 'y ~ x'
Run a clock ~ PCs model for each clock’s DNAmAge and AgeAccel.
var.df <- data.frame(
row.names = c("age", "ageRes", "hannum", "hannumRes", "horvath", "horvathRes", "zhang", "zhangRes", "phenoage", "phenoageRes", "grimage", "grimageRes", "dunedinpaceAge", "dunedinpace"),
variable = c("Age", "AgeAccel", "Hannum", "Hannum residual", "Horvath", "Horvath residual", "Zhang", "Zhang residual", "PhenoAge", "PhenoAge residual", "GrimAge", "GrimAge residual", "DunedinPACE-age", "DunedinPACE"),
var.ex = NaN
)
#Prepare PC labels.
pc.labs <- c("PC1: Neutrophils", "PC2: T cell naive/memory", "PC3: Treg/CD4Tmem", "PC4: Bmem/Mono", "PC5: Bnv/Eos", "PC6: Naive CD8 T cells", "PC7: Mixed (CD4Tnv/CD8Tnv)", "PC8: Eos/Mono", "PC9: Mixed (CD4Tnv+CD8Tmem)", "PC10: Mixed (Bmem)", "PC11: Mixed (Baso)")
PCEffects <- function(form){
fit <- lm(formula = form, data = ss.pc)
var.ex <- round(summary(fit)$r.squared * 100, 1)
print(paste0("Variance explained: ", var.ex, "%"))
#Store variance explained.
var.df[as.character(form)[2], "var.ex"] <<- var.ex
res <- as.data.frame(summary(fit)$coefficients)
colnames(res) <- c("estimate", "std.error", "t.stat", "p.val")
#Calculate 95% confidence intervals.
res$ci.lower <- res$estimate - (1.96 * res$std.error)
res$ci.upper <- res$estimate + (1.96 * res$std.error)
#Select only cell type effects, removing the intercept.
res <- res[(rownames(res) %in% c("PC1", "PC2", "PC3", "PC4", "PC5", "PC6", "PC7", "PC8", "PC9", "PC10", "PC11")),]
res$PC <- factor(rownames(res), levels = rownames(res))
#Add PC labels.
res$PC.label <- pc.labs
res$PC.label <- factor(res$PC.label, levels = unique(res$PC.label))
print(res)
return(res)
}
res.age <- PCEffects(form = formula(age ~ PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC7 + PC8 + PC9 + PC10 + PC11))
## [1] "Variance explained: 43.5%"
## estimate std.error t.stat p.val ci.lower ci.upper
## PC1 0.02332937 0.09791573 0.2382597 8.116918e-01 -0.1685855 0.21524419
## PC2 -6.46156949 0.14944642 -43.2366974 0.000000e+00 -6.7544845 -6.16865452
## PC3 -0.72824707 0.17574455 -4.1437817 3.486108e-05 -1.0727064 -0.38378774
## PC4 0.99617076 0.19500859 5.1083429 3.399769e-07 0.6139539 1.37838760
## PC5 -1.95336043 0.22412537 -8.7154809 4.155524e-18 -2.3926462 -1.51407469
## PC6 -6.28966305 0.22680381 -27.7317344 4.057491e-155 -6.7341985 -5.84512758
## PC7 4.35027432 0.23649668 18.3946525 1.229801e-72 3.8867408 4.81380783
## PC8 -0.44156734 0.24144279 -1.8288695 6.749274e-02 -0.9147952 0.03166052
## PC9 0.83517919 0.26369586 3.1672063 1.550561e-03 0.3183353 1.35202307
## PC10 0.09467097 0.29267373 0.3234693 7.463566e-01 -0.4789695 0.66831147
## PC11 0.76572771 0.31814402 2.4068587 1.613492e-02 0.1421654 1.38928999
## PC PC.label
## PC1 PC1 PC1: Neutrophils
## PC2 PC2 PC2: T cell naive/memory
## PC3 PC3 PC3: Treg/CD4Tmem
## PC4 PC4 PC4: Bmem/Mono
## PC5 PC5 PC5: Bnv/Eos
## PC6 PC6 PC6: Naive CD8 T cells
## PC7 PC7 PC7: Mixed (CD4Tnv/CD8Tnv)
## PC8 PC8 PC8: Eos/Mono
## PC9 PC9 PC9: Mixed (CD4Tnv+CD8Tmem)
## PC10 PC10 PC10: Mixed (Bmem)
## PC11 PC11 PC11: Mixed (Baso)
#Clocks.
res.hannum <- PCEffects(form = formula(hannum ~ PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC7 + PC8 + PC9 + PC10 + PC11))
## [1] "Variance explained: 52.6%"
## estimate std.error t.stat p.val ci.lower ci.upper PC
## PC1 0.6132778 0.08875804 6.909547 5.623812e-12 0.4393120 0.78724355 PC1
## PC2 -7.2567293 0.13546926 -53.567351 0.000000e+00 -7.5222490 -6.99120955 PC2
## PC3 -0.3002849 0.15930783 -1.884935 5.951021e-02 -0.6125283 0.01195843 PC3
## PC4 0.8487059 0.17677018 4.801183 1.634319e-06 0.5022364 1.19517544 PC4
## PC5 -2.0269365 0.20316377 -9.976860 3.557852e-23 -2.4251375 -1.62873548 PC5
## PC6 -6.8969005 0.20559171 -33.546589 6.572806e-218 -7.2998602 -6.49394076 PC6
## PC7 3.6262222 0.21437804 16.915082 4.521235e-62 3.2060412 4.04640313 PC7
## PC8 -0.4476190 0.21886156 -2.045215 4.089820e-02 -0.8765876 -0.01865032 PC8
## PC9 0.6978350 0.23903338 2.919404 3.526367e-03 0.2293296 1.16634047 PC9
## PC10 0.3819960 0.26530106 1.439859 1.499847e-01 -0.1379940 0.90198611 PC10
## PC11 1.2517551 0.28838921 4.340506 1.456290e-05 0.6865122 1.81699795 PC11
## PC.label
## PC1 PC1: Neutrophils
## PC2 PC2: T cell naive/memory
## PC3 PC3: Treg/CD4Tmem
## PC4 PC4: Bmem/Mono
## PC5 PC5: Bnv/Eos
## PC6 PC6: Naive CD8 T cells
## PC7 PC7: Mixed (CD4Tnv/CD8Tnv)
## PC8 PC8: Eos/Mono
## PC9 PC9: Mixed (CD4Tnv+CD8Tmem)
## PC10 PC10: Mixed (Bmem)
## PC11 PC11: Mixed (Baso)
res.horvath <- PCEffects(form = formula(horvath ~ PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC7 + PC8 + PC9 + PC10 + PC11))
## [1] "Variance explained: 43.8%"
## estimate std.error t.stat p.val ci.lower ci.upper PC
## PC1 0.1219737 0.09108413 1.339132 1.806029e-01 -0.05655122 0.3004986 PC1
## PC2 -6.4140141 0.13901952 -46.137507 0.000000e+00 -6.68649232 -6.1415358 PC2
## PC3 -0.4607825 0.16348283 -2.818538 4.847731e-03 -0.78120888 -0.1403562 PC3
## PC4 0.7518009 0.18140282 4.144373 3.477167e-05 0.39625138 1.1073504 PC4
## PC5 -1.5033864 0.20848812 -7.210897 6.601465e-13 -1.91202308 -1.0947497 PC5
## PC6 -5.4686838 0.21097968 -25.920429 3.519637e-137 -5.88220393 -5.0551636 PC6
## PC7 3.3604922 0.21999628 15.275223 2.937058e-51 2.92929951 3.7916849 PC7
## PC8 -0.5768126 0.22459729 -2.568208 1.025815e-02 -1.01702329 -0.1366019 PC8
## PC9 1.0675295 0.24529777 4.351974 1.382481e-05 0.58674585 1.5483131 PC9
## PC10 0.1074964 0.27225384 0.394839 6.929825e-01 -0.42612111 0.6411140 PC10
## PC11 1.0633693 0.29594707 3.593106 3.306486e-04 0.48331303 1.6434256 PC11
## PC.label
## PC1 PC1: Neutrophils
## PC2 PC2: T cell naive/memory
## PC3 PC3: Treg/CD4Tmem
## PC4 PC4: Bmem/Mono
## PC5 PC5: Bnv/Eos
## PC6 PC6: Naive CD8 T cells
## PC7 PC7: Mixed (CD4Tnv/CD8Tnv)
## PC8 PC8: Eos/Mono
## PC9 PC9: Mixed (CD4Tnv+CD8Tmem)
## PC10 PC10: Mixed (Bmem)
## PC11 PC11: Mixed (Baso)
res.zhang <- PCEffects(form = formula(zhang ~ PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC7 + PC8 + PC9 + PC10 + PC11))
## [1] "Variance explained: 46%"
## estimate std.error t.stat p.val ci.lower ci.upper
## PC1 0.10182666 0.09331251 1.0912434 2.752308e-01 -0.08106585 0.28471918
## PC2 -6.54800091 0.14242064 -45.9764891 0.000000e+00 -6.82714536 -6.26885646
## PC3 -0.62004035 0.16748245 -3.7021214 2.166524e-04 -0.94830594 -0.29177475
## PC4 0.90703455 0.18584084 4.8807061 1.097805e-06 0.54278649 1.27128260
## PC5 -1.95795466 0.21358879 -9.1669356 7.560872e-20 -2.37658869 -1.53932064
## PC6 -6.35852571 0.21614130 -29.4183739 1.488056e-172 -6.78216267 -5.93488876
## PC7 4.05406403 0.22537850 17.9878032 1.174807e-69 3.61232218 4.49580588
## PC8 -0.48414463 0.23009208 -2.1041343 3.542828e-02 -0.93512509 -0.03316416
## PC9 0.90632251 0.25129898 3.6065506 3.140400e-04 0.41377650 1.39886851
## PC10 0.07382008 0.27891454 0.2646692 7.912778e-01 -0.47285242 0.62049258
## PC11 0.74454187 0.30318743 2.4557149 1.410223e-02 0.15029451 1.33878923
## PC PC.label
## PC1 PC1 PC1: Neutrophils
## PC2 PC2 PC2: T cell naive/memory
## PC3 PC3 PC3: Treg/CD4Tmem
## PC4 PC4 PC4: Bmem/Mono
## PC5 PC5 PC5: Bnv/Eos
## PC6 PC6 PC6: Naive CD8 T cells
## PC7 PC7 PC7: Mixed (CD4Tnv/CD8Tnv)
## PC8 PC8 PC8: Eos/Mono
## PC9 PC9 PC9: Mixed (CD4Tnv+CD8Tmem)
## PC10 PC10 PC10: Mixed (Bmem)
## PC11 PC11 PC11: Mixed (Baso)
res.phenoage <- PCEffects(form = formula(phenoage ~ PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC7 + PC8 + PC9 + PC10 + PC11))
## [1] "Variance explained: 49.7%"
## estimate std.error t.stat p.val ci.lower ci.upper
## PC1 1.00896092 0.09640441 10.465920 2.600157e-25 0.82000827 1.19791357
## PC2 -7.25540535 0.14713974 -49.309626 0.000000e+00 -7.54379924 -6.96701147
## PC3 -0.24299304 0.17303197 -1.404325 1.602990e-01 -0.58213569 0.09614962
## PC4 0.94273514 0.19199867 4.910113 9.461544e-07 0.56641775 1.31905253
## PC5 -2.53284680 0.22066604 -11.478190 4.929291e-30 -2.96535224 -2.10034137
## PC6 -7.16087308 0.22330313 -32.067947 2.970320e-201 -7.59854723 -6.72319894
## PC7 3.73178703 0.23284640 16.026819 4.280387e-56 3.27540809 4.18816598
## PC8 -0.76464931 0.23771616 -3.216648 1.307192e-03 -1.23057299 -0.29872563
## PC9 0.27279119 0.25962576 1.050709 2.934549e-01 -0.23607530 0.78165768
## PC10 0.08369212 0.28815636 0.290440 7.714946e-01 -0.48109434 0.64847859
## PC11 0.54496334 0.31323353 1.739799 8.197037e-02 -0.06897437 1.15890105
## PC PC.label
## PC1 PC1 PC1: Neutrophils
## PC2 PC2 PC2: T cell naive/memory
## PC3 PC3 PC3: Treg/CD4Tmem
## PC4 PC4 PC4: Bmem/Mono
## PC5 PC5 PC5: Bnv/Eos
## PC6 PC6 PC6: Naive CD8 T cells
## PC7 PC7 PC7: Mixed (CD4Tnv/CD8Tnv)
## PC8 PC8 PC8: Eos/Mono
## PC9 PC9 PC9: Mixed (CD4Tnv+CD8Tmem)
## PC10 PC10 PC10: Mixed (Bmem)
## PC11 PC11 PC11: Mixed (Baso)
res.grimage <- PCEffects(form = formula(grimage ~ PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC7 + PC8 + PC9 + PC10 + PC11))
## [1] "Variance explained: 46.6%"
## estimate std.error t.stat p.val ci.lower ci.upper
## PC1 0.72059896 0.07978605 9.0316407 2.562624e-19 0.56421829 0.8769796
## PC2 -5.57424622 0.12177553 -45.7747638 0.000000e+00 -5.81292626 -5.3355662
## PC3 -0.59797093 0.14320442 -4.1756459 3.033989e-05 -0.87865158 -0.3172903
## PC4 1.08576990 0.15890160 6.8329701 9.561125e-12 0.77432276 1.3972170
## PC5 -2.36538064 0.18262724 -12.9519595 1.269433e-37 -2.72333003 -2.0074312
## PC6 -5.58907745 0.18480975 -30.2423301 2.588131e-181 -5.95130456 -5.2268503
## PC7 2.64687473 0.19270793 13.7351624 5.398072e-42 2.26916718 3.0245823
## PC8 -0.91678279 0.19673823 -4.6599116 3.265219e-06 -1.30238973 -0.5311759
## PC9 0.02666734 0.21487102 0.1241086 9.012355e-01 -0.39447986 0.4478145
## PC10 -0.12448889 0.23848346 -0.5220022 6.016974e-01 -0.59191648 0.3429387
## PC11 0.47650861 0.25923778 1.8381141 6.611887e-02 -0.03159744 0.9846147
## PC PC.label
## PC1 PC1 PC1: Neutrophils
## PC2 PC2 PC2: T cell naive/memory
## PC3 PC3 PC3: Treg/CD4Tmem
## PC4 PC4 PC4: Bmem/Mono
## PC5 PC5 PC5: Bnv/Eos
## PC6 PC6 PC6: Naive CD8 T cells
## PC7 PC7 PC7: Mixed (CD4Tnv/CD8Tnv)
## PC8 PC8 PC8: Eos/Mono
## PC9 PC9 PC9: Mixed (CD4Tnv+CD8Tmem)
## PC10 PC10 PC10: Mixed (Bmem)
## PC11 PC11 PC11: Mixed (Baso)
#Clock residuals (and also DunedinPACE).
res.hannumRes <- PCEffects(form = formula(hannumRes ~ PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC7 + PC8 + PC9 + PC10 + PC11))
## [1] "Variance explained: 21.2%"
## estimate std.error t.stat p.val ci.lower ci.upper
## PC1 0.59118817 0.03324889 17.7806917 3.676070e-68 0.52602036 0.65635599
## PC2 -1.13853297 0.05074697 -22.4354852 3.607497e-105 -1.23799704 -1.03906890
## PC3 0.38926249 0.05967694 6.5228297 7.749891e-11 0.27229570 0.50622929
## PC4 -0.09452752 0.06621836 -1.4275123 1.535094e-01 -0.22431550 0.03526046
## PC5 -0.17737925 0.07610543 -2.3307042 1.981786e-02 -0.32654590 -0.02821260
## PC6 -0.94147537 0.07701494 -12.2245808 8.996195e-34 -1.09242466 -0.79052609
## PC7 -0.49287496 0.08030631 -6.1374371 9.191863e-10 -0.65027534 -0.33547458
## PC8 -0.02951689 0.08198585 -0.3600242 7.188478e-01 -0.19020915 0.13117537
## PC9 -0.09296203 0.08954224 -1.0381919 2.992428e-01 -0.26846481 0.08254076
## PC10 0.29235596 0.09938215 2.9417352 3.282309e-03 0.09756695 0.48714497
## PC11 0.52671879 0.10803100 4.8756262 1.126269e-06 0.31497802 0.73845956
## PC PC.label
## PC1 PC1 PC1: Neutrophils
## PC2 PC2 PC2: T cell naive/memory
## PC3 PC3 PC3: Treg/CD4Tmem
## PC4 PC4 PC4: Bmem/Mono
## PC5 PC5 PC5: Bnv/Eos
## PC6 PC6 PC6: Naive CD8 T cells
## PC7 PC7 PC7: Mixed (CD4Tnv/CD8Tnv)
## PC8 PC8 PC8: Eos/Mono
## PC9 PC9 PC9: Mixed (CD4Tnv+CD8Tmem)
## PC10 PC10 PC10: Mixed (Bmem)
## PC11 PC11 PC11: Mixed (Baso)
res.horvathRes <- PCEffects(form = formula(horvathRes ~ PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC7 + PC8 + PC9 + PC10 + PC11))
## [1] "Variance explained: 5.4%"
## estimate std.error t.stat p.val ci.lower ci.upper
## PC1 0.10132432 0.03720756 2.7232186 6.492667e-03 0.02839750 0.174251133
## PC2 -0.69473001 0.05678900 -12.2335307 8.089432e-34 -0.80603645 -0.583423561
## PC3 0.18380574 0.06678218 2.7523170 5.943903e-03 0.05291266 0.314698818
## PC4 -0.12993279 0.07410244 -1.7534213 7.960543e-02 -0.27517356 0.015307984
## PC5 0.22557796 0.08516669 2.6486642 8.112449e-03 0.05865125 0.392504670
## PC6 0.09844194 0.08618448 1.1422235 2.534287e-01 -0.07047964 0.267363531
## PC7 -0.49003583 0.08986773 -5.4528562 5.252568e-08 -0.66617659 -0.313895071
## PC8 -0.18597116 0.09174723 -2.0269947 4.272840e-02 -0.36579574 -0.006146583
## PC9 0.32829313 0.10020331 3.2762704 1.060748e-03 0.13189465 0.524691603
## PC10 0.02370097 0.11121477 0.2131099 8.312520e-01 -0.19427998 0.241681923
## PC11 0.38560605 0.12089338 3.1896376 1.435391e-03 0.14865504 0.622557070
## PC PC.label
## PC1 PC1 PC1: Neutrophils
## PC2 PC2 PC2: T cell naive/memory
## PC3 PC3 PC3: Treg/CD4Tmem
## PC4 PC4 PC4: Bmem/Mono
## PC5 PC5 PC5: Bnv/Eos
## PC6 PC6 PC6: Naive CD8 T cells
## PC7 PC7 PC7: Mixed (CD4Tnv/CD8Tnv)
## PC8 PC8 PC8: Eos/Mono
## PC9 PC9 PC9: Mixed (CD4Tnv+CD8Tmem)
## PC10 PC10 PC10: Mixed (Bmem)
## PC11 PC11 PC11: Mixed (Baso)
res.zhangRes <- PCEffects(form = formula(zhangRes ~ PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC7 + PC8 + PC9 + PC10 + PC11))
## [1] "Variance explained: 4%"
## estimate std.error t.stat p.val ci.lower ci.upper
## PC1 0.07948475 0.02331971 3.4084791 6.596427e-04 0.033778116 0.125191387
## PC2 -0.35992888 0.03559232 -10.1125445 9.291185e-24 -0.429689814 -0.290167937
## PC3 0.07738236 0.04185551 1.8487976 6.455995e-02 -0.004654432 0.159419158
## PC4 -0.04697150 0.04644345 -1.0113698 3.119000e-01 -0.138000668 0.044057661
## PC5 -0.08727370 0.05337793 -1.6350146 1.021238e-01 -0.191894454 0.017347049
## PC6 -0.33508387 0.05401583 -6.2034379 6.078539e-10 -0.440954906 -0.229212838
## PC7 -0.11207717 0.05632430 -1.9898547 4.667414e-02 -0.222472796 -0.001681545
## PC8 -0.06126741 0.05750227 -1.0654781 2.867232e-01 -0.173971855 0.051437036
## PC9 0.10649377 0.06280208 1.6957044 9.001881e-02 -0.016598315 0.229585850
## PC10 -0.01684377 0.06970348 -0.2416489 8.090645e-01 -0.153462596 0.119775054
## PC11 0.01122495 0.07576951 0.1481461 8.822349e-01 -0.137283294 0.159733203
## PC PC.label
## PC1 PC1 PC1: Neutrophils
## PC2 PC2 PC2: T cell naive/memory
## PC3 PC3 PC3: Treg/CD4Tmem
## PC4 PC4 PC4: Bmem/Mono
## PC5 PC5 PC5: Bnv/Eos
## PC6 PC6 PC6: Naive CD8 T cells
## PC7 PC7 PC7: Mixed (CD4Tnv/CD8Tnv)
## PC8 PC8 PC8: Eos/Mono
## PC9 PC9 PC9: Mixed (CD4Tnv+CD8Tmem)
## PC10 PC10 PC10: Mixed (Bmem)
## PC11 PC11 PC11: Mixed (Baso)
res.phenoageRes <- PCEffects(form = formula(phenoageRes ~ PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC7 + PC8 + PC9 + PC10 + PC11))
## [1] "Variance explained: 19.9%"
## estimate std.error t.stat p.val ci.lower ci.upper
## PC1 0.986125092 0.04216812 23.38556170 1.492941e-113 0.9034756 1.06877460
## PC2 -0.930532010 0.06436018 -14.45819420 3.105188e-46 -1.0566780 -0.80438605
## PC3 0.469847779 0.07568567 6.20788328 5.910684e-10 0.3215039 0.61819168
## PC4 -0.032361369 0.08398186 -0.38533759 7.000076e-01 -0.1969658 0.13224308
## PC5 -0.620810215 0.09652122 -6.43185237 1.407249e-10 -0.8099918 -0.43162863
## PC6 -1.004269455 0.09767470 -10.28177625 1.700205e-24 -1.1957119 -0.81282703
## PC7 -0.526456128 0.10184901 -5.16898618 2.467332e-07 -0.7260802 -0.32683207
## PC8 -0.332423434 0.10397909 -3.19702206 1.399237e-03 -0.5362224 -0.12862442
## PC9 -0.544719574 0.11356253 -4.79664871 1.671518e-06 -0.7673021 -0.32213701
## PC10 -0.008976056 0.12604206 -0.07121476 9.432304e-01 -0.2560185 0.23806638
## PC11 -0.204565200 0.13701103 -1.49305645 1.355004e-01 -0.4731068 0.06397642
## PC PC.label
## PC1 PC1 PC1: Neutrophils
## PC2 PC2 PC2: T cell naive/memory
## PC3 PC3 PC3: Treg/CD4Tmem
## PC4 PC4 PC4: Bmem/Mono
## PC5 PC5 PC5: Bnv/Eos
## PC6 PC6 PC6: Naive CD8 T cells
## PC7 PC7 PC7: Mixed (CD4Tnv/CD8Tnv)
## PC8 PC8 PC8: Eos/Mono
## PC9 PC9 PC9: Mixed (CD4Tnv+CD8Tmem)
## PC10 PC10 PC10: Mixed (Bmem)
## PC11 PC11 PC11: Mixed (Baso)
res.grimageRes <- PCEffects(form = formula(grimageRes ~ PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC7 + PC8 + PC9 + PC10 + PC11))
## [1] "Variance explained: 28.3%"
## estimate std.error t.stat p.val ci.lower ci.upper
## PC1 0.701748692 0.02471451 28.3942038 6.940667e-162 0.65330826 0.75018912
## PC2 -0.353260857 0.03772115 -9.3650595 1.227348e-20 -0.42719432 -0.27932739
## PC3 -0.009543107 0.04435896 -0.2151337 8.296739e-01 -0.09648667 0.07740045
## PC4 0.280858194 0.04922132 5.7060277 1.239179e-08 0.18438441 0.37733197
## PC5 -0.787054164 0.05657056 -13.9127865 5.112909e-43 -0.89793247 -0.67617586
## PC6 -0.506993480 0.05724662 -8.8563048 1.215313e-18 -0.61919685 -0.39479011
## PC7 -0.868171979 0.05969315 -14.5439120 9.473587e-47 -0.98517056 -0.75117340
## PC8 -0.559993834 0.06094158 -9.1890271 6.184599e-20 -0.67943933 -0.44054834
## PC9 -0.648162255 0.06655839 -9.7382508 3.620692e-22 -0.77861669 -0.51770782
## PC10 -0.200983577 0.07387257 -2.7206794 6.542653e-03 -0.34577381 -0.05619334
## PC11 -0.142203784 0.08030142 -1.7708750 7.665676e-02 -0.29959457 0.01518700
## PC PC.label
## PC1 PC1 PC1: Neutrophils
## PC2 PC2 PC2: T cell naive/memory
## PC3 PC3 PC3: Treg/CD4Tmem
## PC4 PC4 PC4: Bmem/Mono
## PC5 PC5 PC5: Bnv/Eos
## PC6 PC6 PC6: Naive CD8 T cells
## PC7 PC7 PC7: Mixed (CD4Tnv/CD8Tnv)
## PC8 PC8 PC8: Eos/Mono
## PC9 PC9 PC9: Mixed (CD4Tnv+CD8Tmem)
## PC10 PC10 PC10: Mixed (Bmem)
## PC11 PC11 PC11: Mixed (Baso)
res.dunedinpace <- PCEffects(form = formula(dunedinpace ~ PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC7 + PC8 + PC9 + PC10 + PC11))
## [1] "Variance explained: 26.3%"
## estimate std.error t.stat p.val ci.lower ci.upper
## PC1 0.92892594 0.04476168 20.75270619 5.446016e-91 0.8411931 1.0166588
## PC2 -1.45073602 0.06831867 -21.23483952 5.796187e-95 -1.5846406 -1.3168314
## PC3 -0.28795801 0.08034073 -3.58420945 3.420878e-04 -0.4454258 -0.1304902
## PC4 0.57841205 0.08914719 6.48828144 9.729261e-11 0.4036836 0.7531405
## PC5 -0.48324924 0.10245778 -4.71656960 2.479421e-06 -0.6840665 -0.2824320
## PC6 -2.11675678 0.10368221 -20.41581422 2.942103e-88 -2.3199739 -1.9135396
## PC7 -0.08865388 0.10811326 -0.82000929 4.122591e-01 -0.3005559 0.1232481
## PC8 -0.86827802 0.11037435 -7.86666512 4.638406e-15 -1.0846117 -0.6519443
## PC9 -0.30841626 0.12054722 -2.55846835 1.054953e-02 -0.5446888 -0.0721437
## PC10 0.01309479 0.13379431 0.09787256 9.220383e-01 -0.2491421 0.2753316
## PC11 -0.14563818 0.14543793 -1.00137690 3.167045e-01 -0.4306965 0.1394202
## PC PC.label
## PC1 PC1 PC1: Neutrophils
## PC2 PC2 PC2: T cell naive/memory
## PC3 PC3 PC3: Treg/CD4Tmem
## PC4 PC4 PC4: Bmem/Mono
## PC5 PC5 PC5: Bnv/Eos
## PC6 PC6 PC6: Naive CD8 T cells
## PC7 PC7 PC7: Mixed (CD4Tnv/CD8Tnv)
## PC8 PC8 PC8: Eos/Mono
## PC9 PC9 PC9: Mixed (CD4Tnv+CD8Tmem)
## PC10 PC10 PC10: Mixed (Bmem)
## PC11 PC11 PC11: Mixed (Baso)
Plot the effect sizes of the clocks versus the effect sizes of the residuals.
plotDNAmAgeAndResidual <- function(dnamage, age.residual, clock.lab){
dat <- rbind(dnamage, age.residual)
dat$type <- factor(c(rep("DNAmAge", 11), rep("AgeAccel", 11)), levels = c("DNAmAge", "AgeAccel"))
dat
#Prepare annotations for variance explained for the clock and its residual.
clock.var <- paste0("Variance explained: ", round(var.df[grep(clock.lab, var.df$variable)[1],"var.ex"], 1), "%")
resid.var <- paste0("Variance explained: ", round(var.df[grep(clock.lab, var.df$variable)[2],"var.ex"], 1), "%")
p <- ggplot(dat, aes(x = PC.label, y = estimate)) +
geom_hline(yintercept = 0, linetype = "dashed") +
geom_point(aes(color = type), shape = 19, size = 1.5) +
geom_errorbar(aes(ymin = ci.lower, ymax = ci.upper, color = type), width = 0.4) +
labs(x = NULL, y = paste0(clock.lab, " effect"), color = "") +
theme_bw() +
scale_color_manual(values = c("#005580", "#993300"), drop = F) +
scale_y_continuous(limits = c(-7.6, 4.7), breaks = seq(-10, 10, 2)) +
theme(axis.text.x = element_text(angle = 45, hjust = 1),
plot.margin = margin(t = 5.5, r = 5.5, b = 5.5, l = 22, unit = "pt"))
print(p)
}
plotDNAmAgeAndResidual(res.hannum, res.hannumRes, "Hannum")
plotDNAmAgeAndResidual(res.horvath, res.horvathRes, "Horvath")
plotDNAmAgeAndResidual(res.zhang, res.zhangRes, "Zhang")
plotDNAmAgeAndResidual(res.phenoage, res.phenoageRes, "PhenoAge")
plotDNAmAgeAndResidual(res.grimage, res.grimageRes, "GrimAge")
#For DunedinPACE, we only have AgeAccel (the clock itself is a putative age acceleration marker). Add a dummy dataframe for the DNAmAge.
res.dunedinpace.dummy <- res.dunedinpace
for(i in 1:7){
if(class(res.dunedinpace.dummy[,i]) == "numeric"){
res.dunedinpace.dummy[,i] <- NaN
} else{
res.dunedinpace.dummy[,i] <- NA
}
}
plotDNAmAgeAndResidual(res.dunedinpace.dummy, res.dunedinpace, "DunedinPACE")
## Warning: Removed 11 rows containing missing values (`geom_point()`).
Make the plots for all clocks in a single faceted figure.
plot.data <- rbind(res.hannum, res.hannumRes, res.horvath, res.horvathRes, res.zhang, res.zhangRes, res.phenoage, res.phenoageRes, res.grimage, res.grimageRes, res.dunedinpace.dummy, res.dunedinpace)
plot.data$type <- factor(rep(c(rep("DNAmAge", 11), rep("AgeAccel", 11)), 6), levels = c("DNAmAge", "AgeAccel"))
plot.data$clock <- factor(c(rep("Hannum", 22), rep("Horvath", 22), rep("Zhang", 22), rep("PhenoAge", 22), rep("GrimAge", 22), rep("DunedinPACE", 22)), levels = c("Horvath", "Hannum", "Zhang", "PhenoAge", "GrimAge", "DunedinPACE"))
head(plot.data)
## estimate std.error t.stat p.val ci.lower ci.upper PC
## PC1 0.6132778 0.08875804 6.909547 5.623812e-12 0.4393120 0.78724355 PC1
## PC2 -7.2567293 0.13546926 -53.567351 0.000000e+00 -7.5222490 -6.99120955 PC2
## PC3 -0.3002849 0.15930783 -1.884935 5.951021e-02 -0.6125283 0.01195843 PC3
## PC4 0.8487059 0.17677018 4.801183 1.634319e-06 0.5022364 1.19517544 PC4
## PC5 -2.0269365 0.20316377 -9.976860 3.557852e-23 -2.4251375 -1.62873548 PC5
## PC6 -6.8969005 0.20559171 -33.546589 6.572806e-218 -7.2998602 -6.49394076 PC6
## PC.label type clock
## PC1 PC1: Neutrophils DNAmAge Hannum
## PC2 PC2: T cell naive/memory DNAmAge Hannum
## PC3 PC3: Treg/CD4Tmem DNAmAge Hannum
## PC4 PC4: Bmem/Mono DNAmAge Hannum
## PC5 PC5: Bnv/Eos DNAmAge Hannum
## PC6 PC6: Naive CD8 T cells DNAmAge Hannum
ggplot(plot.data, aes(x = PC.label, y = estimate)) +
geom_hline(yintercept = 0, linetype = "dashed") +
geom_point(aes(color = type), shape = 19, size = 1.5) +
geom_errorbar(aes(ymin = ci.lower, ymax = ci.upper, color = type), width = 0.4) +
theme_bw() +
scale_color_manual(values = c("#005580", "#993300"), drop = F) +
scale_y_continuous(limits = c(-7.6, 4.7), breaks = seq(-10, 10, 2)) +
theme(axis.text.x = element_text(angle = 45, hjust = 1),
plot.margin = margin(t = 5.5, r = 5.5, b = 5.5, l = 22, unit = "pt")) +
facet_wrap(facets = vars(clock), scales = "free", nrow = 2, dir = "h") +
theme(panel.spacing.x = unit(1.2, "cm")) +
labs(color= "", x = "", y = "")
## Warning: Removed 11 rows containing missing values (`geom_point()`).
Plot variance explained for each clock and residual.
plot.data <- var.df
plot.data$variable <- factor(plot.data$variable, levels = unique(plot.data$variable))
plot.data$group = factor(c("Age", "Age", "Hannum", "Hannum", "Horvath", "Horvath", "Zhang", "Zhang", "PhenoAge", "PhenoAge", "GrimAge", "GrimAge", "DunedinPACE", "DunedinPACE"), levels = c(c("Age", "Hannum", "Horvath", "Zhang", "PhenoAge", "GrimAge", "DunedinPACE")))
plot.data$type = factor(c("Age", "AgeAccel", "DNAmAge", "AgeAccel", "DNAmAge", "AgeAccel", "DNAmAge", "AgeAccel", "DNAmAge", "AgeAccel", "DNAmAge", "AgeAccel", "DNAmAge", "AgeAccel"), levels = c("Age", "DNAmAge", "AgeAccel"))
plot.data$var.round <- paste0(round(plot.data$var.ex, 0), "%")
#Split the data into 1st-gen and 2nd-gen clocks.
plot.data$gen <- c(rep("Age", 2), rep("1st-Generation Clocks", 6), rep("2nd-Generation Clocks", 6))
plot.data$gen <- factor(plot.data$gen, levels = c("Age", "1st-Generation Clocks", "2nd-Generation Clocks"))
#Remove unused rows.
plot.data$var.ex[c(2, 13)] <- 0
plot.data$var.round[c(2, 13)] <- ""
plot.data
## variable var.ex group type var.round
## age Age 43.5 Age Age 44%
## ageRes AgeAccel 0.0 Age AgeAccel
## hannum Hannum 52.6 Hannum DNAmAge 53%
## hannumRes Hannum residual 21.2 Hannum AgeAccel 21%
## horvath Horvath 43.8 Horvath DNAmAge 44%
## horvathRes Horvath residual 5.4 Horvath AgeAccel 5%
## zhang Zhang 46.0 Zhang DNAmAge 46%
## zhangRes Zhang residual 4.0 Zhang AgeAccel 4%
## phenoage PhenoAge 49.7 PhenoAge DNAmAge 50%
## phenoageRes PhenoAge residual 19.9 PhenoAge AgeAccel 20%
## grimage GrimAge 46.6 GrimAge DNAmAge 47%
## grimageRes GrimAge residual 28.3 GrimAge AgeAccel 28%
## dunedinpaceAge DunedinPACE-age 0.0 DunedinPACE DNAmAge
## dunedinpace DunedinPACE 26.3 DunedinPACE AgeAccel 26%
## gen
## age Age
## ageRes Age
## hannum 1st-Generation Clocks
## hannumRes 1st-Generation Clocks
## horvath 1st-Generation Clocks
## horvathRes 1st-Generation Clocks
## zhang 1st-Generation Clocks
## zhangRes 1st-Generation Clocks
## phenoage 2nd-Generation Clocks
## phenoageRes 2nd-Generation Clocks
## grimage 2nd-Generation Clocks
## grimageRes 2nd-Generation Clocks
## dunedinpaceAge 2nd-Generation Clocks
## dunedinpace 2nd-Generation Clocks
# Other visualization order.
plot.data <- plot.data[-c(2, 13),]
ggplot(data = plot.data, aes(x = group, y = var.ex, fill = type, label = var.round)) +
geom_bar(stat = "identity", width = 0.8, position = position_dodge(width = 0.8)) +
geom_text(hjust = 0.5, vjust = 1.15, position = position_dodge(width = 0.8), mapping = aes(color = type), size = 3) +
theme_bw() +
theme(axis.text.x = element_text(angle = 45, hjust = 1), axis.ticks.x = element_blank(), panel.grid.major.x = element_blank()) +
scale_fill_manual(values = c("#222222", "#005580", "#993300")) +
scale_color_manual(values = c("white", "white", "white")) +
scale_y_continuous(breaks = seq(0, 60, 10), limits = c(0, 53)) +
labs(x = "", y = "Variance explained (%)", fill = "") +
guides(color = "none", fill = "none") +
facet_grid(cols = vars(type), scales = "free", space = "free")
Plot variance explained per PC.
PrintTopPCs <- function(){
data.frame(
row.names = paste0("PC", 1:11),
var = rbind(
round(summary(lm(formula = var ~ PC1, data = ss.pc))$r.squared * 100, 3),
round(summary(lm(formula = var ~ PC2, data = ss.pc))$r.squared * 100, 3),
round(summary(lm(formula = var ~ PC3, data = ss.pc))$r.squared * 100, 3),
round(summary(lm(formula = var ~ PC3, data = ss.pc))$r.squared * 100, 3),
round(summary(lm(formula = var ~ PC5, data = ss.pc))$r.squared * 100, 3),
round(summary(lm(formula = var ~ PC6, data = ss.pc))$r.squared * 100, 3),
round(summary(lm(formula = var ~ PC7, data = ss.pc))$r.squared * 100, 3),
round(summary(lm(formula = var ~ PC8, data = ss.pc))$r.squared * 100, 3),
round(summary(lm(formula = var ~ PC9, data = ss.pc))$r.squared * 100, 3),
round(summary(lm(formula = var ~ PC10, data = ss.pc))$r.squared * 100, 3),
round(summary(lm(formula = var ~ PC11, data = ss.pc))$r.squared * 100, 3)
)
)
}
top.pcs <- data.frame(
row.names = paste0("PC", 1:11),
PC = pc.labs,
age = rep(NA, 11),
hannum = rep(NA, 11),
horvath = rep(NA, 11),
zhang = rep(NA, 11),
phenoage = rep(NA, 11),
grimage = rep(NA, 11),
dunedinpaceDummy = rep(NA, 11),
ageRes = rep(NA, 11),
hannumRes = rep(NA, 11),
horvathRes = rep(NA, 11),
zhangRes = rep(NA, 11),
phenoageRes = rep(NA, 11),
grimageRes = rep(NA, 11),
dunedinpace = rep(NA, 11)
)
top.pcs$PC <- factor(top.pcs$PC, levels = (unique(top.pcs$PC)))
#Calendar age
ss.pc$var <- ss.pc$age
top.pcs$age <- PrintTopPCs()$var
#DNAmAge
ss.pc$var <- ss.pc$hannum
top.pcs$hannum <- PrintTopPCs()$var
ss.pc$var <- ss.pc$horvath
top.pcs$horvath <- PrintTopPCs()$var
ss.pc$var <- ss.pc$zhang
top.pcs$zhang <- PrintTopPCs()$var
ss.pc$var <- ss.pc$phenoage
top.pcs$phenoage <- PrintTopPCs()$var
ss.pc$var <- ss.pc$grimage
top.pcs$grimage <- PrintTopPCs()$var
#AgeAccels
ss.pc$var <- ss.pc$hannumRes
top.pcs$hannumRes <- PrintTopPCs()$var
ss.pc$var <- ss.pc$horvathRes
top.pcs$horvathRes <- PrintTopPCs()$var
ss.pc$var <- ss.pc$zhangRes
top.pcs$zhangRes <- PrintTopPCs()$var
ss.pc$var <- ss.pc$phenoageRes
top.pcs$phenoageRes <- PrintTopPCs()$var
ss.pc$var <- ss.pc$grimageRes
top.pcs$grimageRes <- PrintTopPCs()$var
ss.pc$var <- ss.pc$dunedinpace
top.pcs$dunedinpace <- PrintTopPCs()$var
#Remove unused columns.
top.pcs <- top.pcs[,-c(8, 9)]
top.pcs
## PC age hannum horvath zhang phenoage grimage
## PC1 PC1: Neutrophils 0.001 0.559 0.025 0.016 1.362 1.076
## PC2 PC2: T cell naive/memory 26.104 33.622 29.576 28.203 30.225 27.630
## PC3 PC3: Treg/CD4Tmem 0.240 0.042 0.110 0.183 0.025 0.230
## PC4 PC4: Bmem/Mono 0.240 0.042 0.110 0.183 0.025 0.230
## PC5 PC5: Bnv/Eos 1.061 1.166 0.722 1.121 1.638 2.212
## PC6 PC6: Naive CD8 T cells 10.739 13.186 9.335 11.547 12.783 12.061
## PC7 PC7: Mixed (CD4Tnv/CD8Tnv) 4.725 3.353 3.242 4.317 3.193 2.488
## PC8 PC8: Eos/Mono 0.047 0.049 0.092 0.059 0.129 0.286
## PC9 PC9: Mixed (CD4Tnv+CD8Tmem) 0.140 0.100 0.263 0.174 0.014 0.000
## PC10 PC10: Mixed (Bmem) 0.001 0.024 0.002 0.001 0.001 0.004
## PC11 PC11: Mixed (Baso) 0.081 0.221 0.179 0.080 0.038 0.045
## hannumRes horvathRes zhangRes phenoageRes grimageRes dunedinpace
## PC1 6.155 0.173 0.276 10.830 14.280 7.842
## PC2 9.800 3.499 2.427 4.139 1.553 8.211
## PC3 0.828 0.177 0.081 0.763 0.001 0.234
## PC4 0.828 0.177 0.081 0.763 0.001 0.234
## PC5 0.106 0.164 0.063 0.819 3.428 0.405
## PC6 2.910 0.031 0.913 2.093 1.389 7.590
## PC7 0.733 0.695 0.094 0.529 3.747 0.012
## PC8 0.003 0.096 0.027 0.202 1.496 1.127
## PC9 0.021 0.251 0.068 0.456 1.680 0.119
## PC10 0.168 0.001 0.001 0.000 0.131 0.000
## PC11 0.463 0.238 0.001 0.044 0.056 0.018
library(reshape2)
library(ggplot2)
plot.data <- melt(top.pcs, id.vars = "PC")
plot.data$group <- factor(NA, levels = c("Age", "DNAmAge", "AgeAccel"))
plot.data$group[plot.data$variable %in% c("age", "ageRes")] <- "Age"
plot.data$group[plot.data$variable %in% c("hannum", "horvath", "zhang", "phenoage", "grimage", "dunedinpaceDummy")] <- "DNAmAge"
plot.data$group[plot.data$variable %in% c("hannumRes", "horvathRes", "zhangRes", "phenoageRes", "grimageRes", "dunedinpace")] <- "AgeAccel"
#Make a new variable that is the same for DNAmAge and AgeAccel of the same clock for better plotting.
plot.data$clock <- as.character(plot.data$variable)
plot.data$clock[grep("Res", plot.data$clock)] <- do.call("rbind", strsplit(plot.data$clock[grep("Res", plot.data$clock)], split = "Res"))
plot.data$clock[grep("Dummy", plot.data$clock)] <- do.call("rbind", strsplit(plot.data$clock[grep("Dummy", plot.data$clock)], split = "Dummy"))
plot.data$clock <- factor(plot.data$clock, levels = unique(plot.data$clock))
levels(plot.data$clock) <- c("Age", "Hannum", "Horvath", "Zhang", "PhenoAge", "GrimAge", "DunedinPACE")
head(plot.data)
## PC variable value group clock
## 1 PC1: Neutrophils age 0.001 Age Age
## 2 PC2: T cell naive/memory age 26.104 Age Age
## 3 PC3: Treg/CD4Tmem age 0.240 Age Age
## 4 PC4: Bmem/Mono age 0.240 Age Age
## 5 PC5: Bnv/Eos age 1.061 Age Age
## 6 PC6: Naive CD8 T cells age 10.739 Age Age
ggplot(plot.data, aes(x = clock, y = PC, fill = value)) +
geom_tile() +
facet_grid(cols = vars(group), scales = "free", space = "free") +
geom_text(aes(label = format(round(value, 1), 1)), vjust = 0.5, hjust = 0.5, size = 3.2) +
scale_fill_gradient(low = "white", high = "red", limit = c(0, max(plot.data$value)),
guide = guide_colorbar(label.position = "left", label.hjust = 1)) +
theme_bw() +
theme(axis.text.x = element_text(angle = 45, hjust = 1),
axis.title = element_blank(),
panel.grid = element_blank(),
strip.text = element_blank(),
legend.title = element_text(hjust = 1),
legend.position = "left") +
scale_y_discrete(position = "right") +
labs(x = "Clock variable", y = "Cell count PC", fill = "Variance\nexplained\n(%)")
sessionInfo()
## R version 4.3.1 (2023-06-16)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Rocky Linux 8.10 (Green Obsidian)
##
## Matrix products: default
## BLAS/LAPACK: /usr/lib64/libopenblas-r0.3.15.so; LAPACK version 3.9.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## time zone: Europe/Amsterdam
## tzcode source: system (glibc)
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] ggpubr_0.6.0 reshape2_1.4.4 ggplot2_3.4.3
##
## loaded via a namespace (and not attached):
## [1] sass_0.4.7 utf8_1.2.3 generics_0.1.3 tidyr_1.3.0
## [5] rstatix_0.7.2 stringi_1.7.12 lattice_0.21-8 digest_0.6.33
## [9] magrittr_2.0.3 evaluate_0.21 grid_4.3.1 fastmap_1.1.1
## [13] plyr_1.8.8 jsonlite_1.8.7 Matrix_1.6-1 backports_1.4.1
## [17] mgcv_1.8-42 purrr_1.0.2 fansi_1.0.4 scales_1.2.1
## [21] jquerylib_0.1.4 abind_1.4-5 cli_3.6.1 rlang_1.1.1
## [25] munsell_0.5.0 splines_4.3.1 withr_2.5.0 cachem_1.0.8
## [29] yaml_2.3.7 tools_4.3.1 ggsignif_0.6.4 dplyr_1.1.2
## [33] colorspace_2.1-0 broom_1.0.5 vctrs_0.6.3 R6_2.5.1
## [37] lifecycle_1.0.3 stringr_1.5.0 car_3.1-2 pkgconfig_2.0.3
## [41] pillar_1.9.0 bslib_0.5.1 gtable_0.3.4 glue_1.6.2
## [45] Rcpp_1.0.11 highr_0.10 xfun_0.40 tibble_3.2.1
## [49] tidyselect_1.2.0 rstudioapi_0.15.0 knitr_1.43 farver_2.1.1
## [53] htmltools_0.5.6 nlme_3.1-162 labeling_0.4.2 rmarkdown_2.24
## [57] carData_3.0-5 compiler_4.3.1